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Abstract — We present a novel method for computing the 
Minkowski Functionals from isodensity surfaces extracted 
directly from the Delaunay tessellation of a point distribution. 
This is an important step forward compared to the previous 
cosmological studies when the isodensity surface was built 
in the field on a uniform cubic grid and therefore having 
a uniform spatial resolution. The density field representing 
a particular interest in cosmology is the density of galaxies 
which is obtained from the highly nonuniform distribution 
of the galaxy positions. Therefore, the constraints caused 
by the spatially uniform grid put severe limitations on the 
studies of the geometry and shapes of the large-scale objects: 
superclusters and voids of galaxies. Our technique potentially 
is able to eliminate most of these limitations. The method is 
tested with some simple geometric models and an application 
to the density field from an N-body simulation is shown. 

Keywords -DelsiunsLy Tessellations; Minkowski Functionals; 
cosmology 

I. Introduction 

Redshift surveys of galaxies have shown that even on the 
largest observed scales the bulk of matter in the universe is 
concentrated in clusters and superclusters of galaxies which 
are separated by huge almost empty regions, called voids 
of galaxies. Cosmological N-body simulations in the cur- 
rently popular ACDM model reveal that at moderate density 
thresholds J ~ 3, the supercluster network percolates, while 
occupying a small fraction 3%) of the total volume 
(32 1 . At greater density thresholds the structure breaks into 
disconnected parts most of which are clusters of galaxies and 
long filaments that link several clusters of galaxies together 
similar to structures shown in Fig. (H They are often referred 
to as superclusters of galaxies. At even greater thresholds 
only the clusters of galaxies are seen as islands in a sea. 

A number of statistical measures have been suggested 
to quantify the patterns made by spatial concentrations of 
galaxies. The traditional approach to quantify clustering 
makes use of the hierarchy of correlation functions [24J. 
Unfortunately they become cumbersome to evaluate beyond 
the three point function and are not particularly revealing 
as geometry and shapes are concerned. Various statistics 



were proposed to quantify the geometrical and topological 
properties of the structure: e.g. percolation analysis (301, 
1 48 1, counts in cells (TSl . (TVL minimal spanning tree (H, 
the genus measure (T2ll . wavelet analysis (Tsl . Multiscale 
Morphology Filters from the Hessian of the density field (H 
and potential field [14], and recently the SpineWeb algorithm 
(21 and the skeleton of the cosmic web 1211 . l39ll , etc. 

Mecke, Buchert and Wagner (1994) have introduced the 
Minkowski functionals (hereafter MF) to cosmology. The 
four Minkowski functionals computed for the excursion sets 
at many density thresholds contain valuable information 
regarding both the geometrical as well as topological dis- 
tribution of matter in the Universe. There have been three 
major attempts made to study the morphology of the Cosmic 
Web using MF; these efforts differ in their approach of 
evaluating the MF: (1) Firstly, Boolean grain models study 
the MF of surfaces which result due to intersecting spheres 
decorating the input point-set |20|. (2) Secondly, Krofton's 
formulae make it possible to calculate MF on a density field 
defined on a grid |29|. In this case the MF are calculated by 
using the information of the number of vertices, edges, faces 
and cuboids. (3) Finally, an alternative, resolution-dependent 
approach consists in employing the Koenderink invariants 



A radically new and powerful approach for using the 
Minkowski functionals in cosmology was firstly introduced 
for the analysis of two-dimensional CMB (Cosmic Mi- 
crowave Background) maps (22l . (33l . This approach con- 
sists of constructing isotemperature contours and analysis of 
individual parts of the excursion set at many levels. These 
contours are triangulated and the MF are evaluated for the 
resulting closed polygons. Since Minkowski functionals are 
additive in nature one can glean information regarding both 
individual objects (galaxies, clusters, voids) as well as the 
supercluster- void network in its totality. 

However, a three-dimensional algorithm required a new 
component that builds a surface of approximately constant 
(to linear order) density. It was implemented in fTl] and 
then used for the analysis of the mass distribution at the 



non-linear stage in cosmological N-body simulation of the 
ACDM model | 34| and mock galaxy catalogs [351, 136 1. In 
addition, particular ratios of Minkowski functional quan- 
tifies the morphology of large scale structure by telling 
us whether the distribution of matter in superclusters/voids 
is spherical, planar, filamentary etc. |25|. The numerical 
technique used in the cosmological studies is based on 
the marching cubes algorithm (MCA) in three-dimensional 
density field specified on a uniform mesh. It therefore has 
a uniform resolution scale over the region in question. 
The data points are the positions of the galaxies which 
are distributed in highly non-uniform manner: there are 
clusters of galaxies where thousands of galaxies packed in 
the volume of a few Mpc^ and there are voids of galaxies 
where hardly one galaxy can be found in hundreds of Mpc^. 
Thus, for the study of the densest clusters of galaxies one 
needs to have a fine mesh, but it does not work in voids 
were most of sites are empty. On the other hand a mesh 
with large cells would allow to study the voids but erase 
the clusters. A numerical method based on adaptive mesh is 
badly needed for the studies of the structure in the universe. 

A. Extracting implicit surfaces from scalar fields 

The extraction of implicit surfaces from a given scalar 
field is an important problem for visualization and volumet- 
ric data analysis. The marching cubes algorithm |16| was 
introduced in order to extract isosurfaces from a regular grid 
by identifying intersections between each voxel of the grid 
and an isosurface defined by a threshold value 5t and then 
generating a triangular mesh from the intersecting points. 
The MCA is an efficient and fast solution but it presents 
a few ambiguities. In addition to this its use was until 
recently restricted by a patent. In order to circunvent the 
patent restrictions and solve the ambiguities in the MCA 
the Marching Tetrahedra Algorithm (MTA) was developed 
[10], |13|. The main idea behind MTA is to divide each 
voxel of the regular grid into 6 tetrahedra and identify the 
intersections between each tetrahedron and the isosurface 
defined by 5t. This greatly simplifies the identification of 
intersections since there are only 8 cases to consider: no 
intersection, 4 intersections with a triangle and 3 intersec- 
tions with two planar triangles or quad (see Fig. [TJ. There 
are no ambiguities in the intersection cases and they can 
be conveniently encoded in a lookup table. Both MCA 
and MTA use linear interpolation in order to identify the 
intersecting points. 

The MCA and MTA take as input a scalar field sampled 
on a regular grid. In many of the applications of MCA and 
MTA the scalar field is computed directly in a previous step, 
like in tomography and MRI medical images or analytic 
functions for 3D texture and surface generation in computer 
images. In the astronomical context often the scalar field 
(density) must be reconstructed from a point distribution 
such as mass particles from N-body simulations or galaxy 



positions from redshift surveys. The accurate reconstruction 
of such datasets is not trivial an in general one has to rely 
on a series of assumptions about the underlying scalar field 
such as continuity and linearity. 

A major improvement in the reconstruction and analysis 
of density fields from point distributions came with the 
implementation of Voronoi and Delaunay tessellation tech- 
niques |7J, im, ifTn . ll23l . In particular for astronomy 
the Delaunay Tessellation Field Estimator (DTFE) |27J, 
EH, ll44l based on the previous work of O which applied 
the Delaunay tessellation for the interpolation of a field 
sampled by an (in principle) arbitrary point distribution. 
The DTFE estimates local densities at each point as being 
inversely proportional to the volume of the adjacent Voronoi 
cell of the point. The density field is then sampled on a 
regular grid by linearly interpolating the density field inside 
each tetrahedron containing the sampling point. The DTFE 
represents a natural alternative to grid-based methods as it 
adapts itself to the point distribution, it produces a space- 
filling density field and is well defined for all regions inside 
the volume containing the points. However, the DTFE is 
also known to produce strong artifacts as a result of the 
linear interpolation of the density field over very large 
and thin triangles. By construction it can not provide zero 
density estimates, therefore it can not reproduce the density 
field inside very underdense regions where there are no 
sampling points. In its original implementation it samples 
the density field on a regular grid thus the final density field 
is not adaptive and due to the discrete sampling scheme it 
introduces aliasing in the density field. 

The use of a regular grid is not required in the case 
when we already have values of the scalar field at each 
sampling point. By recognizing that the Delaunay interpo- 
lation schema does not require the coupling with the point 
distribution (as in the DTFE) one can use the tessellation 
to interpolate any scalar field independently of the sampling 
points. Here we use a Delaunay tessellation to describe the 
connectivity between points that will be used to interpolate 
the density and extract isosurfaces. The density field can be 
computed by adaptive methods such as the scale-adaptive 
DTFE, Spline or fixed-scale Gaussian kernel. The density 
values at each vertex of the tessellation are then used to 
identify intersections with the isosurface. By computing the 
isosurfaces directly from the Delaunay tessellation we avoid 
two extra steps: i).- interpolation of the density field on a 
regular grid and ii).- tetrahedrisation of the voxels in the 
grid. Our algorithm can be directly applied to the Delaunay 
tessellation with no further processing. 

II. Isosurfaces from the Delaunay Tessellation: 
Marching tetrahedra Algorithm 

The first step in the reconstruction of isosurfaces is 
the computation of the Delaunay tessellation of the point 



distribution which will be used to define the spatial connec- 
tivity between points used to interpolate the density field. 
The tessellation is also used to generate adaptive density 
estimates at each point as described in 1271 

A. Isosurface extraction 

Once we computed the Delaunay tessellation we test each 
tetrahedron for intersection with the isosurface defined by 
the threshold value 5t. We do this by simply evaluating 
mm{5i) < St < max(J^), where (^i=o...3 is the density field 
evaluated at each of the four vertices of the tetrahedron. 
We then proceed to count the number of edge intersections. 
There are 7 intersecting cases (see Fig.[T]), when the isovalue 
divides the tetrahedron in two lower values and two higher 
values the tetrahedron is cut by a quad, otherwise it is cut 
by a single triangle. Each intersecting point corresponds 
to a vertex of the triangulated isosurface. The position of 
the intersecting point is estimated by interpolating the field 
between the two vertices (xq, xi) of the intersected edge as 
follows: 

X = {1 — t) Xq -\- t Xi 

y = {l-t)yo^tyi (1) 
z = {1 -t) Zo^t Zi 

where the parameter t runs in the interval [0, 1] and is given 
by: 

t = {St-So)/{Si-So). (2) 

Depending on the number of intersecting points we construct 
a triangle (3 intersections) or pair of coplanar triangles for a 
quad (4 intersections). The vertices of the triangle are then 
oriented to have their normal vector pointing outside the 
surface i.e. in the direction of decreasing density. This can be 
easily checked with the inner product between the normal of 
the triangle and the less dense vertex of the tetrahedron. The 
orientation of the quads is also checked to avoid incorrect 
triangulations when the vertices are not properly oriented as 
adjacent triangles. 

III. Implementation 

The Delaunay tessellation of the point distribution 
was computed using the publicly available CGAL 
(www.cgal.org) library. This efficient library allows us to 
compute the tessellation for 256^ particles in a couple of 
minutes on a regular linux workstation. Our implementation 
can generate open or periodic boundary conditions by using 
a buffer region around the original box and replicating the 
indexes of each point. Note that the most recent CGAL 
release includes periodic boundaries natively. This new 
feature will be included in the next version of our code. 
For the present application our approach is sufficient. 

The code that constructs the isosurface from the Delaunay 
TesselUation was written as part of a larger visualization 




Figure 1. Possible intersection cases between a tetrahedron and the 
isosurface. a: no intersection with tetrahedron, b-d: intersection in a quad, 
there are two more dense vertices and two less dense, e-h: intersection in 
a triangle, there is one single denser or less dense vertex. 

suite written in C++ built on top of the OpenFrameworks 
(http://www.openframeworks.cc) environment for OpenGL- 
based volume and surface renderings as well as multimedia 
mixing for data sonification, sound visualization, etc. (see 
Fig. El). 

A version of the MTA was also written in the IDL 
language for fast development and testing of the algorithms. 
The code loops over all tetrahedra and constructs triangles or 
quads as described in section III This raw surface consists of 
unconnected triangles so we proceed to remove duplicated 
vertices by indexing them. Repeated vertices are identified 
by means of the edge from which they were created. Al 
vertices sharing the same parent edge correspond to the same 
point. 

A. Isolated Object Extraction 

In general, the isosurface consists of a set of isolated 
objects. In order to study the properties of each individual 
object we must be able to identify it. This is done by first 
constructing a list of adjacent triangles. For each triangle 
in the isosurface we iteratively link and label its three 
adjacent triangles until all triangles are connected. The final 
triangulation divided into isolated objects is stored in a 
binary format and as a Waveform .obj format file that can 
be used for analysis and visualization. 

B. Density estimation 

In order to interpolate the density field inside each tetra- 
hedron we need density estimates at each point. There are 
several methods for this purpose, among them the n-nearest 
neighbor, top-hat counts, Gaussian kernel and DTFE are the 
most common (O, ll38l , B2l . Adaptive density estimators 
such as the DTFE are strongly dependent on the spatial 
sampling. This is more evident in low density regions where 
one can see strong artifacts arising from the linear interpola- 
tion inside the tetrahedron. The use of averaging procedures 
on the tessellation or higher order interpolation alleviates 



some artifacts at the expense of increased complexity of the 
method. The end result is that isosurfaces extracted from 
adaptive density estimates tend to be highly irregular in low 
density regions. For the analysis presented in the coming 
sections we choose a fix-scale Gaussian kernel density 
estimator. This limits out ability to probe into compact high 
density regions such as the interior of clusters. However, 
for low-intermediate density regions the Gaussian smoothing 
produces better surfaces than the adaptive approaches like 
the DTFE provided that the smoothing kernel is larger 
than the mean interparticle separation. By using a Gaussian 
kernel we still take advantage of the space-filling properties 
of the tessellation and the convenient use of tetrahedra 
to interpolate the density field and identify intersections 
with the isosurface. In a following paper we will present a 
more detailed comparison between adaptive and fixed- scale 
density estimation. For the purposes of this work we restrict 
our analysis to Gaussian densities. 

We compute Gaussian densities at particle i as the sum 
of a Gaussian kernel evaluated at each of the surrounding j 
particles: 



(3) 



where a is scale of the smoothing kernel. For convenience 
we use the overdensity defined as: 



P 



(4) 



Where p is the mean density computed for the given a. 

The density estimation is an 0{N'^) operation, which 
makes the brute-force implementation unfeasible for large 
datasets. We use a bucket algorithm where the containing 
box is divided into a coarse grid. Instead of evaluating the 
Gaussian kernel for all the particles we restrict the operation 
to particles inside bucket containing particle i as well as its 
26 adjacent buckets. We use a grid size such that there is 
a warranty of enclosing particles up to 5cr away from the 
target particle. Our fully threaded code greatly improves the 
speed of the density estimation and scales almost linearly 
with the number of processors. 

The isodensity surfaces computed assuming a linear vari- 
ation of the density field inside each tetrahedron tend to be 
very irregular. This effect is even more pronounced in the 
interface regions between high and low density cosmological 
structures such as the edges of walls, filaments and clusters 
where the density field profile is close to a power law. In 
these regions, if the spatial sampling is low, there can be 
a large variation in the density field inside an individual 
tetrahedron. Also, the overall density distribution of the 
evolved Cosmic Web is roughly log-normal 1 8 1 . As a result 
of this, the isodensity surfaces computed from linear inter- 
polation tend to be highly irregular and to overestimate the 
enclosed volume inside the surface. Based on the character 



of the density distribution and the power-law behavior of the 
density profiles it makes more sense to use a transformation 
that linearizes the field inside the tetrahedron. We therefore 
use the log{S) and linearly interpolate the log-transformed 
density field. The log {6) interpolation does not conserve 
mass as in the linear case but it produces smoother surfaces 
and better volume estimation. 

IV. Morphological Parameters 
A. Minkowski Functionals 

In this work we discuss the geometry and topology of 
the regions bounded by the isodensity surfaces and therefore 
make no prior assumptions about the shapes of superclusters 
and voids. It is worth noting that some regions may have 
more than one boundary surface and possess nontrivial 
topology of the boundaries. The complete characterization of 
an arbitrarily complex region in three dimensions obviously 
cannot be achieved if only a few numbers are used. At best 
one can try to design some basic characteristics that serve a 
particular purpose. Our purpose is to provide basic measures 
suitable for quantification of typical components of the large- 
scale structure: superclusters and voids. 

Four Minkowski functionals are effective non-parametric 
descriptors of the morphological properties of surfaces in 
three dimensions |[T9ll , (201, E3- They are 

• Volume V enclosed by the surface S, 

• Area A of the surface, 

• Integrated mean curvature C of the surface. 
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da, 



(5) 



where Ri and R2 are the principal radii of curvature 
at a given point on the surface. 
the Euler characteristic 
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(6) 



B. Shapefinders 

As demonstrated in 1251 , ||26l particular ratios of 
Minkowski functionals called shapefinders provide us with 
a set of non-parametric measures of sizes and shapes of 
objects. Therefore, in addition to determining MFs we shall 
also derive the shapefinders, T (Thickness), B (Breadth) and 
L (Length) defined as follows: 

This normalization gives roughly a half of the diameter 
in the case of anisotropic objects. The three shapefinders 
describing an individual region bounded by one or several 
isolated surfaces of constant density have dimensions of 
length and provide us with an estimate albeit quite crude of 
the region's 'extensions': T is the shortest and thus describes 
the characteristic thickness of the region or object, L is 



T 



(7) 



typically the longest and characterizes the length of the 
object; B is intermediate and can be associated with the 
breadth of the object. This simple interpretation is obviously 
relevant only for fairly simple shapes. The choice of the 
coefficients in (|7]) results in a sphere having all three sizes 
equal to its radius T = B = L = R. A triaxial ellipsoid 
has values of T, B and L close but not equal to the lengths 
of its three principal semi-axes: shortest, intermediate and 
the longest respectively. It is worth noting that T, B and L 
are only the estimates of three basic sizes (semi-axes) of an 
object which work quite well on such objects as a triaxial 
ellipsoid and torus 1251 , ll26ll , 1371 but no three numbers can 
describe an arbitrary, complex three-dimensional shape. In 
particular, the length of a region with many tunnels may 
be better characterized by L = — x/2)), because 

(1 — x/2) is equal to the number of tunnels and therefore 
L approximately represents the length of the parts between 
tunnels. 

As suggested in 1251 a rough idea about the global 
'shape' of a region can be provided by an additional pair 
of dimensionless shapefinders: 

where P and F are measures of Planarity and Filamentarity 
respectively (P, F < 1). A sphere has P = F = 0, an 
ideal filament has P = 0, P = 1 while P = 1, P = for an 
ideal pancake. Other interesting shapes include 'ribbons' for 
which P P 1. When combined with the genus measure, 
the triplet {P, P, G} provides an example of shape-space 
which incorporates information about topology as well as 
morphology of superclusters and voids. 

Note however that non-geometrical shapefinding statistics 
based on mass moments etc. can give misleading results 
when applied to large scale structure, as demonstrated in 

Ei. 

C Numerical Estimates of Minkowski Functionals 

We construct a 2-polytopal surface using an assembly 
of triangles in which every triangle shares its sides with 
each of its three neighboring triangles. Here we outline the 
technique used for the numerical estimate of the Minkowski 
Functionals. It was suggested in (371 . where it is described 
in more detail, see also |36l 

• The area of such a triangulated surface is 

Nf 

A = Y,A,, (9) 

where Ai is the area of the i^^ triangle face and Np 
is the total number of faces which compose a given 
surface. 



• The volume enclosed by this 2-polytopal surface is the 
summed contribution from Np tetrahedra 

Nf 

y = 

i=l 

Vi = ^MnjP^)i. (10) 

Here Vi is the volume of an individual tetrahedron 
whose base is a triangle on the surface. (rijP^) is the 
scalar product between the outward pointing normal n 
to this triangle and the mean position vector of the three 
triangle vertices, for which the j^^ component is given 

-. 1 . . . 

p^ = -{Pl + pi + pi). (11) 

The subscript i in ([TOl) refers to the i—th tetrahedron, 
while the vectors Pi , P2 , P3 in (fTTl) define the location 
of each of three triangle vertices defining the base of 
the tetrahedron relative to an origin. The origin can be 
arbitrary chosen. 

• The extrinsic curvature of a triangulated surface is 
localized in the triangle edges. As a result the integrated 
mean curvature C is determined by the formula 

C= -e, (12) 

where iij is the length of the edge common to triangles 
i and j and is the angle between the normals fii & 
fij to these adjacent triangles 

cos^^j = rii ■ rij. (13) 

The summation in ([T2b is carried out over all edges. It 
should be noted that for a completely general surface, 
the extrinsic curvature can be positive at some (convex) 
points and negative at some other (concave) points on 
the surface. 

• The Euler characteristics of a polyhedral surface is 
given by 

X = Nf-Ne^Nv, (14) 

where Np, Ne^ Ny are, respectively, the total number 
of faces, edges, and vertices defining the surface. For- 
mula [141 is correct for arbitrary shaped faces, obviously 
it can be reduced to 

X = Nv- Nf/2, (15) 

in the case when all faces are triangles since the number 
of edges is Ne = 3Nf/2, which also holds true for 
polytopal surfaces. 

V. Results 

In this section we present two set of tests focused on 
simple geometries and a more complex surface extracted 
from a cosmological N-body simulation. 




Figure 2. Four simple geometries used for testing, from left to right: sphere, pancake, filament and torus. 



A. Simple Models 

We constructed a set of simple geometries in order to test 
our method. Three geometries represent the generic types of 
the gravitational collapse in the Zel'dovich framework (STJ, 
BtI : a compact clump represented by a sphere, pancake and 
filament. Additionally we construct a torus in order to test its 
topology. The geometries were constructed by generating a 
set of hypothetical sampling points over a continuous scalar 
field defined as: 



d. 



1/2 



(16) 



Where d is the scalar value of the field and we define the 
following limiting cases: 



a = b = c 
a ^b^ c 
a^b ^ c 



sphere 

pancake 

filament. 



(17) 



The torus was computed as: 



d={^[R-{x'^y'f^')\z'^ 



1/2 



(18) 



Where R is the radius of the torus and t is the radius of 
the tube. We evaluated the distance function for each of the 
sampling points and proceed to construct the isosurface at 
a defined threshold value. Since the isosurfaces defined on 
the distance field enclose underdense regions we flip the 
normals so that they point outside the surface. Note that in 
this test surfaces the sampling points are not coupled with 
the scalar field. They are used only to define the tessellation 
that will be used to interpolate the density field. 

Fig. |2] shows the four models after we extraced the 
isosurface. We rendered the raw surface (i.e. no normal 
smoothing) in order to emphasize the triangles that form the 
surface. The MTA produces very smooth surfaces although 
in this case it only reflects the spatial sampling of the random 
points. 



Table U sumarizes the result of our analysis over the four 
geometries. The parameters for which there are analytical 
formulae {e.g. all parameters for the sphere) were computed 
with the accuracy better than 1%. The numbers of faces, 
edges and vertices in the case of the sphere were 8,553, 
25,653 and 17,102. The numbers of the triangulation param- 
eters in the cases of the filament, pancake and torus were 
similar while two largest superclusters in the cosmological 
N-body simulation they were about and order of magnitude 
greater. 

Table I 

Minkowski functionals and shape measures 



Mink Func 


Sphere 


Pancake 


Filament 


Torus 


Volume 


0.065 


0.016 


0.004 


0.059 


Area 


0.788 


0.449 


0.163 


1.183 


C 


3.156 


2.564 


1.826 


5.923 


EC 


2 


2 


2 





Thickness 


0.248 


0.109 


0.075 


0.149 


Breadth 


0.249 


0.175 


0.089 


0.199 


Lenght 


0.251 


0.204 


0.145 


0.471 


Planarity 


0.002 


0.233 


0.086 


0.143 


Filamentarity 


0.002 


0.075 


0.238 


0.404 



B. Cosmological Density Field 

We applied our method to a more realistic and complex 
dataset consistinng of a dark-matter N-body simulation of 
a standard ACDM cosmology. The simulation is contained 
inside a periodic box of 200 Mpc. Initial conditions 
were generated on a 256^ grid with Vtm = 0.3, = 0.7, 
(78 = 0.9 and = 73 km s~^ Mpc~^.0 After having set up 
the initial conditions, we follow the subsequent gravitational 
evolution to the present time using the public N-body code 
Gadget2 [|40 I . The simulation box is large enough to contain 
the full range of environments in the Cosmic Web, from 
large underdense voids to massive dense clusters. 

^One megaparsec (Mpc) corresponds to 3.26 x 10^ light years, Qrn and 
Qa correspond to the ratio of matter and dark energy relatively to a close 
universe with ^totai = 1- The Hubble constant H measures the expansion 
of the Universe as a recession velocity per Mpc. 



Figure 3. Volume rendering of a density field (left) and isodensity surface at 6 = 3 (right). 



Table II 

Minkowski functionals and shape measures of two largest 
superclusters in the n-body simulation 



Mink Func 


Green 


Red 


Volume 


0.00359 


0.00202 


Area 


0.429 


0.234 


C 


13.34 


6.54 


EC 


-10 


-6 


Thickness 


0.0251 


0.0258 


Breadth 


0.0322 


0.0358 


Lenght 


1.06 


0.521 


Planarity 


0.124 


0.162 


Filamentarity 


0.941 


0.871 



From the final particle distribution we estimate the local 
density using a Gaussian kernel of 2 Mpc/h. We then 
compute the Delaunay tessellation with periodic conditions 
and applied the MTA (section |ll|) in order to construct an 
isodensity surface corresponding to an overdensity threshold 
of (5 = 3. From the raw isosurface we label individual objects 
as described in section IIII-AI and compute their MF and 
shapefinders (see section lIV-Ab . Fig. [3] shows the volume 
rendering of the density field inside the simulation box (left 
panel) and the isosurface extracted at S = 3 (right panel). 
The isosurface is formed by a large number of unconnected 
structures enclosing overdense regions and clumps of matter. 

Fig. 131 shows the two largest superclusters depicted from 
the N-body simulation. The largest structure is shown in 
blue and the second largest in red. Their parameters, MF, 
thickness, breadth, length and shapefinders are given in Table 
HH Both have comparable thicknesses and breadths while the 
length of the first largest supercluster about twice longer than 
that of the second largest cluster. The superclusters have 6 
and 4 tunnels respectively, a couple of which can be seen 
in the figure. Although the figure gives some idea about the 



complexity of the objects, its full appreciation is possible 
only with the help of visualization software where the three- 
dimensional structure of each can be fully recognized. 

VI. Conclusion and future work 

We present a novel technique for computing Minkowski 
functional i.e. the volume, area, integrated mean curvature 
and Euler characteristic of complex objects (superclusters 
and voids of galaxies) that form due to the nonlinear 
gravitational growth of initial Gaussian random fluctuations 
in dark matter density field. These structures can be reli- 
ably modeled in cosmological N-body simulations and are 
observed in large galaxy redshift catalogs like SDSS and 
2dF. The technique can also be used for other applications 
when quantitative information is needed for arbitrary shaped 
structures. An important step forward has been made by 
replacing the Marching Cube Algorithm in previous studies 
1 34 1, 1 35 1, 1 37 1 by Marching Tetrahedra Algorithm for 
building the isodensity surface. The MCA operates on the 
fields specified on a uniform cubic grid and therefore has 
a uniform spatial resolution. The MTA finds the density 
field directly from the positions of galaxies or particles in 
cosmological N-body simulations by means of Delaunay 
tessellation and therefore has an adaptive spatial resolution. 
Potentially the galaxy density field can be analyzed with 
better resolution in the regions where more galaxies are 
observed or more simulation particles are clustered however 
this requires additional testing which under way and will be 
reported in the following papers. 

In this work we report the results of the basic tests that 
shows that the code works and ready for the studies of the 
density fields obtained in cosmological N-body simulations 
of the structure in the universe or real structures in galaxy 
redshift catalogs. 



Figure 4. Largest (blue) and second largest (red) structures at 6t = 3 and Gaussian smoothing of 2 Mpc/h. 
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